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1 .0 SUMMARY 

In this document, a method is developed for using the integrals of systems of 
nonlinear, ordinary, differential equati<ms in a numerical integration pro- 
cess to (1) control the local errors in these integrals and (2) reduce the 
global errors of the solution. The method is general and can be applied to 
either scaijjr or vector integrals. A number of example problems , with 
accompanying numerical results, are used to verify the analysis and support the 
conjecture of global error reduction. 


2.0 INTRODUCTION 

Whittaker (ref. 1) defines the integral of a system of differential equations as 
a function of the state and time when the total time derivative is zero and the 
state variables are any functions of time. In solving nonlinear differential 
equations by a numerical integration process, these integrals will be used along 
with their associated numerical errors to (1) control the local errors in these 
integrals and (2) reduce the global error of the solution. 

Integrals of systems of differential equations are constraints on the solution 
and, as such, can be used to reduce the number of degrees of freedom in the prob- 
lem. Topologically, as more integrals of a system are introduced into the prob- 
lem, the solution is constrained to lie on a larger manifold of the solution 
space with a resulting reduction in the global errors. The limiting case is, of 
course, an analytic solution to the differential equations where the solution is 
known at any time and the integral and global error are zero. 

The direct approach (analytical substitution) to using integrals of systems of 
differential equations to reduce the number of degrees of freedom in a prob- 
lem and to reduce the errors introduces ether types of difficulties (such as sin- 
gularities and switching logic in the remaining unsolved equations). Invariably, 
the overhead of calculating the right-hand side (RHS) of the remaining differen- 
tial equations increases significantly; thus, any overall advan age of such an ap- 
proach is nullified. 

In a direct analytical approach, Szebehely (ref. 2) linearized certain nonlinear 
differential equations by utilizing the integrals of the system and by intro- 
ducing a new independent variable. This was an effort to formalize the results 
of Steifel and Scheifele (ref. 3), Burdet (ref. 4), Szebehely (ref. 5), and 
Sperling (ref. 6) who attempted to linearize, and in some cases stabilize, cer- 
tain nonlinear differential equations by transformations of the dependent and in- 
dependent variables of the problem. 

A direct numerical approach v«?as made by Nacozy (ref. 7). The numerical errors 
in the integrals of the system were used to rectify the solution at each integra- 
tion step in an attempt to stabilize the solution and reduce the errors . For 
some integrals, a linear expansion was necessary to compute the correction 
vector. Using a fourth-order predictor -corrector integration routine with a vari- 
able stepsize, global error reduction of two or three orders of magnitude v^as 
obtained . 
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An optimization technique is the moat general, indirect approach to controlling 
the error in a system described by nonlinear differential equations. A perfor- 
mance function is defined as a function of the error in the Integrals, and the 
differential equations of the system are adjoined to this performance function 
by Lagrange multipliers. Through the use of variational calculus or some other 
similar technique, differential equations for the multipliers and optimality and 
boundary conditions can be developed. However, this is not an advisable pr'oce- 
dure for controlling the errors for the following reasons: (1) the number of de- 
grees of freedom in the problem typically increases twofold, (2) an iteration 
procedure is introduced to satisfy the developed boundary conditions, a’ld (3) 
the overhead in the solution of the problem increases significantly. However, 
the error in the integral would be minimized and would be indirectly used to re- 
duce the global error of the solution. 

A functional, indirect approach to introducing Integrals of systems in the solu- 
tion of differential equations was advanced by Baumgarte. In a number of 
studies (refs. 8 throu^ 10) he used the Integrals of the system to stabilize 
certain nonlinear differential equations and reduce the errors. The procedure 
was based on the principle of adjoining a form of the constraint (the coeffi- 
cient of the second derivative) by a Lagrange multiplier to the original sec- 
ond-order system of differential equations. The technique was applied to sev- 
eral problems and the results were encouraging. The error in the integral 
(almost always the energy) or constraint was substantially reduced by using this 
control process. However, several characteristics of these studies were 
disappointing: (1) global error results of the solution were almost always 

absent although some analytic solutions were available, (2) the technique re- 
quired a partlculeu? formulation for the problem, (3) the parameters that were 
introduced were not mathematically defined, and (4) a lack of generality existed 
when applying the approach to any system of equations and constraints. 

This document does not pretend to advance the best technique for using and 
controlling the errors in the Integrals of a system of differential equations. 

It does propose a general technique for incorporating the integrals of a system 
of nonlinear differential equations and their associated errors in a process 
that will control the integral errors and reduce the global error of the solu- 
tion. The process requires no increase in the number of degrees of freedom in 
the solution. The technique has two disadvantages: (1) it requires some 

premathematical analysis to formulate the control vector, and (2) it generates 
some additional overhead and complexity in the solution. 


3.0 STATEMENT €F PROBLEM 

This study examines the numerical solution of a first-order, nonlinear system 
of ordinary differential equations of the form 


X = F(X,t) X(0) = Xq at t = 0 


( 1 ) 


which possess integrals of motion of the form 
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J(X,t) s K (2) 


where X and F are n-vectors and J and K are m-vectors; m < n and K 
is a constant defined by the Initial conditions. The system of equation (1) 
is assumed not to be analytically integrable in terms of known functions 
and, hence, must be solved by a numerical integration process. Specifically, 
this method presents a numerical solution to equation (1) with the additional 
criteria that the solution lies "arbitrarily close" to the (m + 1) -dimensional 
manifold described in equation (2). 

One obvious solution to this problem is to use 'Equation (2) to eliminate m of 
the X's and thus reduce the problem to integrating only the remaining (n-m)- 
first-order differential equations. However, past experiences and the examina- 
tion of a particular system of differential equations and their associated inte- 
i^als indicate that this is not an advisable procedure. There is information 
about the solution embedded in equation (2), and it should not be used only as 
a check on the numerical integration of equation (1), but it should be used 
either directly or indirectly in the solution to reduce the errors and possibly 
the number of degrees of freedom in the problem. For example, if a problem- free 
procedure could be devised for introducing the m Integrals of motion into a 
system of n differential equations, and the m Integrals were introduced into 
the solution one at a time, the number of degrees cf freedom in the solution 
would be reduced by one each time, requiring the solution to lie on a manifold 
with a dimension increasing by one. If there were actually n integrals of mo- 
tion, then as m approaches n, the solution would be constrained to a larger 
subspace of the problem, and would also have global errors that tend to zero- 
vanish when m equals n. 

Because the direct approach to using the Integrals of motion in the solution of 
a system of differential equations introduces other mathematical difficulties, 
an indirect approach shall be proposed to solve equation (1) while attempting to 
satisfy the constraints expressed by the integrals of the motion. 


4.0 SOLUTION WITH ERRORS 

If equation (1) is solved by a numerical integration process, the solution will 
certainly contain errors. Consider these errors as due to the inability of the 
numerical integration process to correctly evaluate the right-hand side (RHS) of 
the differential equations. Defining the numerical errors in the RHS from the 
integration process as 5F, the differential equation (eq. (1)) can be expressed 
as 


X = F(X,t) + 6F (3) 


where at the initial time (t = 0), 6F(0) = 0 . The vector F(X,t) in equation 

(3) is the exact representation of the RHS of the differential equations. Now 
consider the term 6F as a perturbation to the original system of equation (1). 
The integrals of motion (eq. (2)) are not conserved (K is not equal to a 
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constant) and the differential equations for their rate of change can be 
expressed as 


K = G(X, 6F, t) 


(^) 


where G is a vector function of the indicated arguments. The numerical 
error in the integral of motion is defined as 


e = K - Ko , Kq = K( 0 ) 


(5) 


The differential equation for the time rate of change of the error is simply 



( 6 ) 


However, 6F is not known (except at the initial time) since the exact solution 
of the RHS in equation (1), F(X,t) is not known. Hence, it appears that this 
development is an Interesting but Insignificant exercise. 


5.0 SOLUTION WITH CONTROL 

In this section, a solution philosophy used in optimal control theory (ref. 11) 
will be adopted. In controlling a system, it is fundamental to have a process 
that is (1) observable and (2) controllable. The first criterion is cer- 
tainly fullfilled, for it is noted in the numerical integration process that the 
value of the integral is not constant but grows in some manner characteristic to 
the particular numerical integrator, stepslze, etc. The second criterion, how- 
ever, is not fullfilled. 

The process is not controllable because the error vector 6F is an unknown 
output of the numerical integration process and not an input. 

A control vector X (an n-vector) is added to the RHS of aquation (1) (the 
exact equation) in an attempt to control the numerical error fiF of the integra- 
tion process. 

X = F(X,t) + X (7) 


The justification for the addition of the vector X to the original equation 
(eq. (1)) is: if m integrals of the motion or constraints exist to a system 
of n equations, then the solution to this system of equations is constrained 
to lie on a m-dimensional manifold of the solution space and there are only 
n-m independent degrees of freedom in the problem. The introduction of 
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a control vector X is an attempt to numerically reduce the n-order d'?pendent 
aystem of equations to a (n-m)-order independent system of equations. 

the introduction of the control vector X to the RHS of the exact equation (eq. 
0)) implies that the numerically integrated equation (eq. (3)) has a similar 
term added to its RHS. Of course, it is understood that with the addition of 
the control vector X, the term 6F appearing in equation (3) will have a dif- 
ferent meaning since the error vector will certainly be disturbed by the intro- 
duction of this control vector. The introduction of the control vector X is 
an attempt to introduce a term to cancel all or part of the error vector 6F such 
that when the numerical integration process is applied to equation (7), the 
error (defined in eq , (5)) will be arbitrarily close to zero. Thus, if 


6F + X s 0 (8) 


and equations (5) and (6) are used to obtain a "stable" solution for the 
control vector X, the numerical integration process will produce a value 
of the state X, which nulls the error. Also, since the number of degrees 
of freedom in the system of equations to be Integrated have been numerically 
reduced, it is conjectured that the global error of the solution will also be 
reduced. A solution for the n-vector X from the system of m-constralnt 
equations must now be obtained. 

Since the error rate and the error are now controllable, a stable differential 
equation for the desired functional relationship between these two errors 
is Introduced as 


e = -Ye ( 9 ) 


where Y is a positive function (defined in the appendix). Now, any error 
arising in the integral of motion due to a dissatisfaction of equation ( 5 ) 
will be critically damped by the control vector X obtained from equation 
(9). Using equations (M), (5), (6), and (8) in equation (9) gives 


G(X,-X,t) = -Y (K(X,t) - Kq) (10) 


Since the vector X must span-the-space defined by the state vector X, a 
solution for X is assumed to be 
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where A is an undefined matrix of appropriate dimensions. Using equation (11) 
in equation (10) yields 


G(X, ^X,t) = -Y (K(X,t) Ko) 


( 12 ) 


The problem of controlling the integral errors is reduced to determining a solu- 
tion of an algebraic equation. Any solution for the undefined matrix A that 
satisfies equation (12) will produce a control vector (from eq. (11)) that, 
when used in the numerical integration process, will control the errors in 
the integrals of motion. The difficulty in determining a matrix 
A that satisfies equation (12) is, of course, problem dependent. In two of the 
example problems (linear oscillator and two-body problem) with a scalar Integral 
of motion (energy), the matrix A was obtained by inspection. In a third exam- 
ple (two-body problem with an angular momentum integral), some manipulation was 
required to obtain a matrix A that satisfied equation (12). In a final exam- 
ple, a solution is developed to the two-body problem when both the energy and 
the angular momentum Integral errors are present. 


6.0 ONE-DIMENSIONAL HARMONIC OSCILLATOR 

The first-order linear, differential equations describing the one-dimensional 
harmonic oscillator state are 

Xi = X2 

X2 = -Xi (13) 


Since an analytic solution to this problem exists , there are two integrals of mo- 
tion. For this exercise, however, it is assumed that equation (13) cannot be 
solved analytically and that only one integral of motion (ene^rgy) exists and is 
defined as 


J(X) = 1/2 X^x = k 


where the superscript refers to the transpose. 
The numerical error in the integral is defined as 


e = 1/2 X^X - kc 


ko = k(0) 


(14) 
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Adding the control vector X to equation (13), the eofstpolled equation to 
be integrated is 


Xi = X2 + Xi 

X2 = -Xi + X 2 (15) 


Developing the total time derivative of the integral of notion, using equa- 
tion (15) yields 


k = G(X,X) = X^X ( 16 ) 


A stable differential equation for the functional relationship between 
the error and error rates is defined as 


e = -Y e (17) 


where Y is a positive scalar function. From equations (14), (16), and (17), 
the control vector X is required to satisfy the following equation. 


X^(X + - X) = Y ko 
2 ° 


( 18 ) 


If Y were a constant, equation (18) would represent a new integral of motion; 
however, one that is functionally dependent on the control vector, 

A solution is assumed for the control of the form 


X + - X = Yko aX (19) 

2 


where a in an undefined coefficient. Using equation (19) in equation (I 8 ), 
a necessary condition is 


a X^x = 1 


f 
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One value of the coefficient satisfying the equation is 


Qt 




which results in a control vector (from eq. (19)) of 


X 


Y e X 
i k 


(20) 


As expected, the control vector for the energy integral error is in the form 
of a feedback control law directly proportional to the error. It should 
be noted that the vector X is never singular unless the energy constant 
is zero (trivial case). 


6.1 COMPARATIVE RESULTS 

Equation (15) was integrated with a fixed step, fourth order, Runge-Kutta 
integrator using the control vector defined by equation (?0). Since the period 
of the uncontrolled solution is equal to 2 tt, solutions were obtained for a con- 
stant integration stepsize defined as 


h = 2 tt/N 


where N is the number of steps in the integration process . The function Y 
was determined from a solution of the equation e(t + h) = 0 at each integration 
step (see appendix) . For this linear problem, the function Y for each value 
of N was found to be a constant. 

In the solution, it was noted that the uncontrolled error in the energy grew 
in a linear manner with time. The controlled error remained essentially 
constant at a value five to ten orders of magnitude less (depending on the 
value of N and t) than the uncontrolled error. Thus, the method described 
in the analysis of controlling the error in the integral of motion appears 
to be valid for the linear problem. 

Because the analytic solution to equation (13) is known, the global error of the 
numerically integrated solution can also be computed. The global error at a 
given time is defined as 


AX = Ix^ - X^l 


( 21 ) 
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where the superscripts I and A on the vectov X refer to the Integrated and 
analytic values, respectively. The global error of the controlled solution 
ms found to be almys less than that of the uncontrolled solution. This indi- 
cates that the correct information from the error in the integral of motion 
ms entering the solution via the control vector X. However, although the 
integral errors were reduced substantially by the control, the global error 
showed only an infinitesimal reduction. This agrees with the results reported 
by Nacosy (ref. 7) for the harmonic oscillator problem. 

Conclusion: For this linear and stable problem, controlling the error in the in" 

tegral of motion has only a negligible effect on the global error of the solu- 
tion . 


7.0 TWO-BODY PROBLEM WITH ENERGY INTEGRAL 


The first-order, nonlinear differeutial equations describing the two-body prob- 
lem are 


R = V 



(22a) 

(22b) 


where R and V are the position and velocity vectors (respectively), P is 
the gravitational constant and r = | r| . There are three integrals of motion to 
this system of equations (not all independent): two vector integrals - Laplace 
and angular momentum, and one scalar integral - energy. Ilils example is 
concerned only with the energy integral that is formally obtained by scalar 

multiplying equation (22a) by V, and equation (22b) by R, then taking the 
difference and noting the exact differential. This integral can be expressed 
as 


J(R,V) = 1/2 V^V - y/r = k , (23) 


where k is a constant of the motion defined by the initial conditions. If 
equatlcx: (22) is solved with a numerical integrator, the solution will contain 
errors. This numerical error in the energy integral is defined as 


e = k - ko , ko = k(0) ( 2 ^) 


and the differential equation for its time rate of change is defined as 
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• • 
e s k 


Now a control vector 


is added to the RHS of equation (22) in an attempt to control the numerical 
error in the integral of motion. For convenience of notation, a state 
vector is defined as 



and an associated matrix is defined as 



where 0 and I are appropriate null and Identity matrices, respectively. 
Using these identities, the numerical error and error rate can be expressed as 


e = sTaS - ko 

e = -X'^S (25) 


The stable differential equation (eq. (17)) is now used for the functional 
relationship between the error in the integral of motion and its time rate 
of change. Using equation (25) in equation (17) yields 


ST(X - Y AS) = -Yko (26) 

Since the vector X must span the space defined by the state vector S, 
an assumed solution for X is 


X - Y AS = -Yko A S 


(27) 
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Where A is an undefined matrix. Using equation (27) in equation (26), 
a necessary condition is 


sTAs = 1 


(28) 


Conparing equations (28) and (2^), it is noted that one solution for the undefined 
matrix A is 


A 



which results in a control vector (from eq. (27)) of 



(29) 


Thus, the control vector for the energy Integral of the two-body problem can be 
oast in the same form as that of the control vector for the harmonic oscillator 
problem. It is only singular for the special case when the energy constant is 
zero (parabolic orbit). 


7.1 NUMERICAL RESULTS 

Equation (22) with the added control vector defined by equation (29) was in- 
tegrated with the same processor as the previous example. However, the constant 
stepslze was determined as 



where P is the period of the orbit defined by the Initial conditions. 

Figure 1 illustrates a comparison of controlled and uncontrolled solutions for 
a two-body problem. The solutions were obtained using a fourth-order, Runge- 
Kutta integrator with a constant stepsize of N. The abscissa label of NORBIT 
refers to the number of orbits the equations were integrated, while the ordinate 
is the logarithm to the base 10 of either the absolute value of the position 
error or velocity error. (These two errors were approximately equal for this 
problem.) Hence, tlie ordinate is a measure of the number of signifioant digits 
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of Information remaining In the solution at a given time (a measure of the 
global error of the solution) . The solid curves in figure 1 refer to solutions 
obtained with the addition of the control vector X, while the dashed curves 
refer to no control. 

In figure 1(a), the initial conditions were specified by a circular orbit that 
1s an orbit with an eocentrioity, e = 0 and semimajor axis a = 1. For N = 20, 
the uncontrolled solution has lost all information content at NORBIT equals 20, 
whereas the controlled solution still retains approximately two digits of ac- 
curacy . For N s 1<0 , the controlled solution has approximately two aditional 
digits of accuracy at NORBIT equals 40. For a given N, the two solutions are 
diverging; that is, the uncontrolled solutions have a steeper slope or are 
approaching zero more rapidly than the controlled solutions. 

In figure 1(b) and 1(c), similar results are obtained for initial orbits with 
eccentricities of e = 0.1 and e = 0.2, respectively . For N = 20, the 
uncontrolled solution for e = 0.1 and e = 0.2 has zero-significant digits of 
accuracy at NORBIT = 15 and 9, respectively. 

When overlaying figures 1(a) and 1(b), it is noted that the controlled solution 
shows almost no loss in accuracy for this change in eccentricity, whereas the 
uncontrolled solution shows significant losses. Similar (but not as 
significant) results are obtained when overlaying figures Kb) and 1(c). 

The. errors in the integral of motion (energy) in figuras 1(a) through 1(c) 
behave similar to that of the harmonic oscillator problem. The uncontrolled in- 
tegral error grows almost linearly with time; the slope is determined by the ini- 
tial conditions and the number of steps N. The controlled integral error re- 
mains nearly constant in value, with a number of orders in magnitude less than 
the uncontrolled solution. The actual number of orders in magnitude depends on 
the initial conditions, integration stepsize, and integration time. 

In figure 2, the Y function is plotted versus time for one orbit for the 
three sets of initial conditions defined by e = 0.0, 0.1, and 0.2. For e = 0.0, 
which is similar to the previous linear problem, the Y function is nearly 
constant. The Y functions for e ^ 0 are periodic in nature. For e = 0.2, 
there is a portion of the solution where Y < 0. Since there are no constraints 
imposed on either the sign or the magnitude of the function Y, a negative 
Value is certainly possible (although disconcerting) when viewed from a solution 
to equation (9). 

In the next two example problems ( concerned in part with controlling the error in 
the angular momentum vector), only the control vector and an outline of the 
solution are developed. 


8.0 TWO-BODY PROBLEM WITH ANGULAR MOMENTUM INTEGRAL 

The differential equations to be integrated are the same as in the previous exam- 
ple; however, the angular momentum integral shall now be considered, and 
it is expressed as 
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J (R,V) = R X V s K 


(30) 


where the vector K is a constant of the motion defined by the initial con- 
ditions. Similar to the previous example, the numerical error vector 
in this Integral is defined as 


e = K - Ko , Ko = K(0) (31) 

T T 

and after introducing a control vector = (Xr , Xy) into the RHS of equation 
(22), the time rate of change of the error may be determined as 


c=RxXv-VxXr 


(32) 


A stable vector differential equation for the functional relationship between 
the vector error and its rate is defined as 



(33) 


where y is again a positive function. The control vector X must now span a 
space defined by the vectors R, V, and K or 


'or ay ok' 


■r' 

X = 


V 

.3r 3r 3k. 


,K. 


(34) 


where the a's and 
any component of the 
any component of the 


3's are undetermined scalars. But 
vector Xr parallel to the vector 
vector Xy parallel to the vector 


from equation (32), 
V and , similarly , 

R will have no 


effect on the error rate vector e. 

Hence, the vectors Xr and Xy have the simpler form 


Xr = cxR R + OK K 

Xy = 3y V + 0K K (35) 

Using equation (35) in equation (32), the error rate vector is 

e = (ttR + 3y) K + (3{f R - aK V) x K ( 36 ) 
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The first and second terms on the right-hand side of this equation are the rates 
of change of the angular momentum error vector along and perpendicular to the 
vector K, respectively. Taking the Inner product of the vectors K, U s R x K 
and W B V X K (respectively) with equation (36) and using equation (33) 
yields the necessary conditions on the coefficients o and 0: 


(or + By) = -"'^q 

Bk uTu - Or uTw = -Y U^e (37) 

Bk - Or W^W s -Y W^e 




4 


Where 


eT K 

ic^ 


The first equation Implies that oir and By are homogenous in the term q. 

One solution for these coefficients is 

( ) 

«R = -Y c q , By = -Y(1 - c)q (38 


where the coefficient c is to be determined. The latter two equations (37) 
can be solved for Ur and Br as 


“K = r [(e'^w)(uTu) - (e'*’u)(uTw)] 

8k = - [(e'^W)(uTw) - (eTu)(wTw)] 

A 

where the determinant, A, is 
A = (uTu)(wTw) - (uTw)2 


(39) 


It should be noted that the determinant is only zero for rectilinear motion. 
The control vector for the vector angular momentum error is 


^ ) 
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Xr = ap R + Ok K 

Xy = 3v V + $K 


(40) 


where o and B are defined by equations (38) and (39) and Y is determined from 
the condition |e(t h| = 0 (see appendix). 

The coefficient c in the angular momentum error control vector is a scaling pa- 
rameter between the components of error along and normal to the vector K 
(eq. (38)). Since the physics of the problem dictates that the error normal to 
the vector K will be extremely small, the actual value of the coefficient c 
used in the solution will have a minimal effect on the global error. One 
difficulty should be mentioned. Since the angular momentum error is formed 
by a vector product rather than as in the energy error (a scalar product) a con- 
siderable amount of algebraic manipulation will be required to obtain the 
function Y. 


9.0 TWO-BODY PROBLEM WITH ANGULAR MOMENTUM AND ENERGY INTEGRALS 

In this section the scalar energy integral is added to the vector angular momen- 
tum integral, and a control vector X for both of these errors is determined. 
From equations (28) and (37), the necessary conditions on the coefficients a 
and . 8 are 


+ Bv v'^V = Y e. 
r k k 


otR + 8y = -Yjj q 


(41) 


vdiere aK 3 k given by equation (39) and the subscripts k, K refer 

to variables associated with the energy and angular momentum Integrals, respec- 
tively. Solving equation (41) for the coefficients aR and Bv yields 


«R = -(Yj^ ^k ^K VTy)/(2k + p/r) 

e ^ 0 

Bv = (Yj^ ^/r)/(2k + p/r) 


(42a) 


However, this solution is singular vAien the eccentricity of the orbit 
e = 0. A solutioi valid for e = 0 is 
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«R 


. Tk 

k 


0V 


2k 


e s 0 


(i*2b) 


with— the-addJLt.tx)nal. condition 


Yk ek 
2k 


* Q 


(43) 


Thus, for orbits with e = 0, the function Yj{ is not directly determined by 
the integrator but by the function and the integral errors. But, the co- 

efficients given by equation (42b) will produce the energy error control vector 
given by equation (29 )• Then if the coefficients ai( and 3 k are small and 
if e ~0, the control vector that was used to reduce the energy integral error 
will also reduce the angular momentum integral error. I^ie conclusion has been 
numerically verified. 

FoP the problem illustrated in figure 1(a) with N = 20, the energy integral 
control Vector reduced the angulaur momentum integral errors by approximately 
five orders of magnitude. However, when the eccentricity e was equal to 0.1, 
the energy Integral control vector only reduced the angular momentum Integral 
errors by approximately one order of magnitude. 


10.0 RECOMMENDATIONS 

This method for using the Integrals of systems of nonlinear differential equa- 
tion and their associated numerical errors should be applied to (1) a variable 
step integrator, preferably the Runge-Kutta 4/5, (2) other problems where inte- 
grals or other type of constraints are satisfied through rectification at each 
integration step, and (3) an unstable system of nonlinear differential 
equations. 
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APPENDIX 

CALCULATION OF THE FUNCTION Y 

For a fourth-order, fixed-step Runge-Kutta integrator, the state vector is 
updated at a time t + h by the expressim 

X(t + h) = Xo + ? (f1 + 2(F2 + f3) + F**) , Xo*X(t) (A1) 

o 

where h is the integration step and the functions fJ are 
defined as 

f1 s F(Xq, to) 

f2 = ♦ 5 F>. t, * 5 ) 

f3 = f(Xo * 5 f 2, to * 5) 

F** = F(Xq + hF3, to + h) 

For the controlled 'solution, the right-hand sides of the differential equations 
may be separated into two parts; the vector F and the vector r; adjoined by 
the unknown function Y ; 

X s F(X,t) + Y n(X,t) 


where y = X* For simplicity of notation, the argugirxt t is excluded 
from the vectors F and q- Then the vector F^ may be expressed as 


f 1 = F(Xo) + Yn(Xo) = Fo + Y ho U2) 

and the vector F^ may be expressed as 

f2 = F (Xo + ^(Fo + Y no)) + Y n (Xo + ^(Fo + y ho)) 


But the vector \ represents a small perturbation to the vector F, and hence the 
vector f 2 may be approximated as 
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p2 = Fi + Y (ni + 6i) + ... (A3) 

where 

Fl = F(Xi) , m = n(Xl) , Xi = Xo + ^ Fo 
end 

* - - Y IJPAI + fio) * 6o = 0 

^ ^ ax/ xi 

Similarly, the vectors P3 and F^ may be expressed as 
p3 = F 2 + Y(r )2 6?.) 

, (AH) 

F** = F 3 + Y(n 3 + 63 ) 

where 

X2 = Xo + ^ Fi 

X 3 = Xo + h F 2 


Thus, the coefficient of 


h 

6 


in equation (A1) may be expressed as 


pi + 2 (f2 + f3) + F^ = y(7t+,/) 


(A6) 


where 


= Fq + 2 (F-| + F2) + F3 
%' = no + 2(ni + n2) + ns 
s 60 + 2(61 + 62) 63 


For clarity, the analysis is restricted to a particular problem - the harmonic 
oscillator. Extensions to other problems are straightforward. The desire 
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Is to determine the function Y, which will result at a time t h In an 
Integral error of zero. 

0 = e(t + h) 5 1/2 X(t + h)T x(t + h) - ko (A6) 

Using equations (A1) and (A5) In equation (A6) yields 

1 Xo + ? [*5^- (2 = 2ko (A7) 

6 

the function Y; however, If the small term 
(A7) may be expressed as the quadratic equation 

0 

where 

b = |(x * Sj-) ’■ 

o = I X . - 2ko 

There are two solutions for the function Yj however, It may be readily verified 
that the plus sign of the radical Is correct. 


Tl^ls Is a poly nominal In 
3n 

~ Is Ignored, equation 
dX 

a Y2 + b Y + c = 


NASA-JSC 
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